*goal				Data Analysis: Management of Epilepsy in Southwest Rural China								*
*author				Hongmei Yi, Email: hmyi.ccap@pku.edu.cn											*
*date				April 13, 2020																*
*input data			epilepsy.dta													*
*do	file			analysis.do	
*Update				by Hongmei Yi, August 8, 2020


capture		clear
clear		matrix
clear		mata

set			more off
query		memory
set			maxvar	30000	
set			matsize 400


global		user	""		// working directory 

cd			"$user"
use			"$user\epilepsy.dta",clear

//define  heuristic decision-maker
gen 		heuristic = (rot_numofd == 1)
label		var heuristic "Heuristic decision-maker, 1=yes"

count

global		clinician1 male minority age college physician train_epilepsy heardofcpw heuristic
global		clinician1_c age 
global		clinician1_b male minority college physician train_epilepsy heardofcpw heuristic

global		clinic1 numofdoc monthlyvisit InEqu_rich dis_town dum_numofeplp_18
global		clinic1_c numofdoc monthlyvisit dis_town 
global		clinic1_b InEqu_rich dum_numofeplp_18

global		process vig440_r_qe_num  vig440_r_qe_pct vig440_e_qe_num vig440_e_qe_pct vig440_r_q_*  vig440_e_q_* vig440_r_e_* vig440_e_e_* 
global		diagnosis vig440_d_correct vig440_d_pcorrect vig440_d_incorrect
global		advice vig440_yz_regularvisit vig440_yz_drugstaken vig440_yz_life vig440_yz_mental vig440_yz_parents vig440_yz_exercise vig440_yz_safety
global		management 	vig440_corrmng ///
						vig440_pcorrmng vig440_corrmngA1 vig440_corrmngB1 ///
						vig440_incorrmng ///
						vig440_drugpres   ///
						vig440_corronly_cond vig440_unneonly_cond vig440_harmonly_cond  ///
						vig440_corrunne_cond  vig440_corrharm_cond vig440_unneharm_cond ///
						vig440_corrharmunne_cond  
						
			
order	form_code $clinician1 $clinic1 $process $diagnosis $management $advice

tab 	nationality

*Table 1 
tabstat	*certified,stats(count sum mean)

tabstat		$clinician1_c $clinic1_c if physician == 0,stats(median iqr) format(%9.0f)
ci		proportions $clinician1_b $clinic1_b if physician == 0,wilson //% and ci
tabstat	$clinician1_b $clinic1_b if physician == 0, stats(sum)	//n

tabstat		$clinician1_c $clinic1_c if physician == 1,stats(median iqr) format(%9.0f)
ci		proportions $clinician1_b $clinic1_b if physician == 1,wilson
tabstat	$clinician1_b $clinic1_b if physician == 1, stats(sum)	//n

tabstat		$clinician1_c $clinic1_c,stats(median iqr) format(%9.0f)
ci		proportions $clinician1_b $clinic1_b ,wilson
tabstat	$clinician1_b $clinic1_b , stats(sum)	//n

foreach var of varlist		$clinician1 $clinic1 {
	ttest	`var', by(physician) 
	}

ranksum	college,by(physician)


*Table 2
tabstat		$process, stats(median iqr) format(%9.2f)
tabstat		$process if physician == 1, stats(median iqr) format(%9.2f)
tabstat		$process if physician == 0, stats(median iqr) format(%9.2f)

foreach var of varlist		$process {
	ttest	`var', by(physician) 
	}
ttest		vig440_r_qe_pct,by(physician)
ttest		vig440_e_qe_pct,by(physician)
	
*Table 3
ci		proportions $diagnosis ,wilson
tabstat	$diagnosis, stats(sum)	//n

ci		proportions $diagnosis if physician == 1 ,wilson
tabstat	$diagnosis if physician == 1, stats(sum)	//n

ci		proportions $diagnosis if physician == 0,wilson
tabstat	$diagnosis if physician == 0, stats(sum)	//n

foreach var of varlist		$diagnosis  {
	ranksum	`var', by(physician) 
	}

*Table 4
ci		proportions $management ,wilson
tabstat	$management, stats(sum)	//n

ci		proportions $management if physician == 1 ,wilson
tabstat	$management if physician == 1, stats(sum)	//n

ci		proportions $management if physician == 0,wilson
tabstat	$management if physician == 0, stats(sum)	//n

foreach var of varlist		$management  {
	ttest	`var', by(physician) 
	}

gen  vig440_corronly_uncond= vig440_corronly_cond
replace	vig440_corronly_uncond = 0 if vig440_corronly_uncond == .
ci	proportions vig440_corronly_uncond, wilson
	
*Table 5
ci		proportions $advice ,wilson
tabstat	$advice, stats(sum)	//n

ci		proportions $advice if physician == 1 ,wilson
tabstat	$advice if physician == 1, stats(sum)	//n

ci		proportions $advice if physician == 0,wilson
tabstat	$advice if physician == 0, stats(sum)	//n

foreach var of varlist		 $advice  {
	ttest	`var', by(physician) 
	}
					
*supplementary Table 1
asdoc	sum	vig440_r_q1-vig440_r_q15,  stats(mean sd sum count) replace ///
						title(Supplementary Table ST1 Checklist of the Recommended Questions) ///
						save($user\tableA1.doc) fhc(\b) ///
						fs(12) font(Times Newman) dec(4)
*supplementary Table 2						
asdoc	sum	vig440_r_e1-vig440_r_e6,  stats(mean sd sum count) replace ///
						title(Supplementary Table ST2 Checklist of Clinical Examinations and Laboratory Tests)) ///
						save($user\tableA2.doc) fhc(\b) ///
						fs(12) font(Times Newman) dec(4)					

		
*Table 6

kdensity 	monthlyvisit,normal
sum			monthlyvisit,detail
replace		monthlyvisit = 0.1 if monthlyvisit == 0
gen 		ln_monthlyvisit = ln(monthlyvisit)
label		var ln_monthlyvisit "Log form of outpatient visits in the previous month, person"
kdensity	ln_monthlyvisit,normal
sum			ln_monthlyvisit,detail //偏度和标准差在1:1左右比较符合正态分布


gen 		cid = int(form_code/1000)
tab			cid,gen(cid)

global		cfix cid2 cid3 cid4 cid5 cid6 cid7 cid8 cid9 cid10
global		clinic2 numofdoc ln_monthlyvisit InEqu_rich dis_town dum_numofeplp_18
global		clinician2 male minority age college  physician train_epilepsy heardofcpw heuristic 
 
global		p_adjusted $clinician2 $clinic2 $cfix

*define outcomes
gen 		diag = (vig440_diagnosis == 1|vig440_diagnosis == 2)
replace		diag = . if missing(vig440_diagnosis) == 1
label		var diag "Correct or partially correct diagnosis, 1=yes"
gen 		managed = (vig440_corrmng == 1|vig440_pcorrmng == 1)
replace		managed = . if missing(vig440_corrmng) == 1|missing(vig440_pcorrmng)==1
label		var managed "Correct or partially correct management, 1=yes"

reg			vig440_r_qe_pct $p_adjusted, cluster(cid)
outreg2		using "$user\table6.doc", word replace label  ///
			 alpha(0.01, 0.05, 0.1) ///
			 title(Table 6 Determinants of Correct or Partially Correct Diagnosis and Management of Juvenile Epilepsy) ///
			 ctitle(Process quality) ///
			 bdec(2) cdec(3) ///
			 drop($cfix) ci  ///
			 sortvar(heuristic )
			 
logit		diag vig440_r_qe_pct $p_adjusted, cluster(cid) 
mfx,at(mean)
outreg2		using "$user\table6.doc", word append label  ///
			 alpha(0.01, 0.05, 0.1) ///
			 title(Table 6 Determinants of Correct or Partially Correct Diagnosis and Management of Juvenile Epilepsy) ///
			 ctitle(Diagnosis) ///
			 bdec(2) cdec(3)  ///
			 drop($cfix) mfx ci ///
			 sortvar(heuristic vig440_r_qe_num)

logit		managed	diag vig440_r_qe_pct $p_adjusted, cluster(cid)
mfx,at(mean)
outreg2		using "$user\table6.doc", word append label  ///
			 alpha(0.01, 0.05, 0.1) ///
			 title(Table 6 Determinants of Correct or Partially Correct Diagnosis and Management of Juvenile Epilepsy) ///
			 ctitle(Management) ///
			 bdec(2) cdec(3) mfx ci ///
			 sortvar(heuristic vig440_r_qe_num diag)

*Supplementary Table ST5: Effect of RCT
*Appendix effect of treatment
gen 		itt = (grouptype == "A")
reg			vig440_r_qe_pct itt 
logit		diag itt 
mfx,at(mean)  
logit		managed	itt  
mfx,at(mean)  


**Figure 2
gen 			r540_diag = (r540_diagnosis == 1| r540_diagnosis == 2)
label	var 	r540_diag "Inclusion of correct or partially correct diagnosis made by heuristic process, 1=correct, 2=p correct"

tabstat		diag,by(r540_diag) stats(mean n sum)
gen			vig440_r_qe_num_d = vig440_r_qe_num
replace		vig440_r_qe_num_d = 5 if vig440_r_qe_num>5
replace		vig440_r_qe_num_d = 1 if vig440_r_qe_num<=1

tabstat		diag if vig440_r_qe_num_d== 1,by(r540_diag) stats(mean n sum)
tabstat		diag if vig440_r_qe_num_d== 2,by(r540_diag) stats(mean n sum)
tabstat		diag if vig440_r_qe_num_d== 3,by(r540_diag) stats(mean n sum)
tabstat		diag if vig440_r_qe_num_d== 4,by(r540_diag) stats(mean n sum)
tabstat		diag if vig440_r_qe_num_d== 5,by(r540_diag) stats(mean n sum)

ranksum		diag if vig440_r_qe_num_d == 1,by(r540_diag)
ranksum		diag if vig440_r_qe_num_d == 2,by(r540_diag)
ranksum		diag if vig440_r_qe_num_d == 3,by(r540_diag)
ranksum		diag if vig440_r_qe_num_d == 4,by(r540_diag)
ranksum		diag if vig440_r_qe_num_d == 5,by(r540_diag)

clear

